!
!  match    finding of cross identification between two lists
!
!  Copyright © 2010-2012 F.Hroch (hroch@physics.muni.cz)
!
!  This file is part of Munipack.
!
!  Munipack is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!  
!  Munipack is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!  
!  You should have received a copy of the GNU General Public License
!  along with Munipack.  If not, see <http://www.gnu.org/licenses/>.
!
!
!  http://en.wikipedia.org/wiki/Backtracking
!
!  derived codes by N.Wirth: Data Structures + Algorithms = Programs
!

module matcher

  implicit none

  integer, parameter, private :: rp = selected_real_kind(15)
  integer, private, dimension(:), allocatable :: id1, id2, fid1, fid2
  real(rp), private, dimension(:), pointer :: x1,y1,x2,y2
  real(rp), private :: pmin, fp0, dmin
  integer, private :: n0, fn0
  integer, private :: nmatch = 7

  type item
     real(rp) :: p
     integer, dimension(:), allocatable :: id1,id2
  end type item
  integer, parameter, private :: maxitems = 100
  integer, private :: nitems
  type(item), dimension(:), allocatable :: items

contains

subroutine match(x_1,y_1,x_2,y_2,nm,id_1,id_2,p_0)

  implicit none

  real(rp), dimension(:), intent(in), target :: x_1,y_1,x_2,y_2
  integer, dimension(:), allocatable, intent(out) :: id_1, id_2
  integer, intent(in) :: nm
  real(rp), intent(out) :: p_0
  integer :: i,i1,i2,j1,j2
  logical :: q
  character(len=80) :: fmt
  real(rp) :: a,amin,dc,dcmin

  x1 => x_1
  y1 => y_1
  x2 => x_2
  y2 => y_2

  nmatch = nm

  if( size(x1) /= size(y1) .or. size(x2) /= size(y2) ) then
     write(*,*) 'Array dimensions does not matches.'
     stop 41
  end if

  ! improper set of counts of objects in sequences
  if( .not. (nmatch > 0) ) stop 42 !.and. maxmatch >= nmatch) ) stop 42

  allocate(id1(nmatch),id2(nmatch))
  allocate(fid1(1),fid2(1))
  nitems = 0
  allocate(items(maxitems))
  do i = 1, maxitems
     allocate(items(i)%id1(nmatch),items(i)%id2(nmatch))
  end do

  fn0 = 0
  fp0 = 0.0
  dmin = huge(dmin)

  do i1 = 1, size(x1) - 2
     do i2 = i1 + 1, size(x1) - 1

        id1(1:2) = (/i1, i2/)

        do j1 = 1, size(x2)
           do j2 = 1, size(x2)

              if( j1 /= j2 ) then
                 id2(1:2) = (/j1, j2/)
                 n0 = 2
                 pmin = 1.0_rp
!                 write(*,*) id1(1:2),'->',id2(1:2)
                 call quest(i2+1,q)
              end if

           end do
        end do

     end do
  end do

  if( fn0 > 0 ) then

     allocate(id_1(fn0),id_2(fn0))
     id_1(1:fn0) = fid1(1:fn0)
     id_2(1:fn0) = fid2(1:fn0)
     p_0 = fp0

  else

     allocate(id_1(1),id_2(1))
     p_0 = 0.0

  end if

  amin = huge(amin)
  dcmin = huge(dcmin)
  write(fmt,"(a,i0,a,i0,a)") "(f7.3,a,",fn0,"i3,a,",fn0,"i3,f10.4,f10.6)"
  do i = 1, nitems
     a = angels(x1,y1,items(i)%id1,x2,y2,items(i)%id2)
     dc = dist(x1,y1,items(i)%id1,x2,y2,items(i)%id2)
     write(*,fmt) items(i)%p,':',items(i)%id1(1:fn0),' ->',items(i)%id2(1:fn0),a,dc
     if( a < amin .or. dcmin > dc) then
!     if( dcmin > dc ) then
!        write(*,*) a,dc
        amin = a
        dcmin = dc
        id_1(1:fn0) = items(i)%id1(1:fn0)
        id_2(1:fn0) = items(i)%id2(1:fn0)
        p_0 = items(i)%p
     end if
  end do

  write(*,*) p_0,id_1,id_2

  do i = 1, maxitems
     deallocate(items(i)%id1,items(i)%id2)
  end do
  deallocate(items)

  deallocate(id1,id2,fid1,fid2)

end subroutine match

recursive subroutine quest(i3,q)

  integer, intent(in) :: i3
  logical, intent(out) :: q
  integer :: i,i1,i2,j1,j2,j3,imax
  real(rp) :: u1, v1, u2, v2, puv, pmax
  character(len=80) :: fmt

  q = .false.

  if( locate(i3,id1(1:n0)) ) return

  i1 = id1(n0-1)
  i2 = id1(n0)
!  write(*,*) i1,i2,i3
  call triangle(x1(i1),y1(i1),x1(i2),y1(i2),x1(i3),y1(i3),u1,v1)

  do j3 = 1, size(x2)

     q = .false.

     if( .not. locate(j3,id2(1:n0)) .and. n0 < nmatch .and. id1(n0) /= i3 ) then
        j1 = id2(n0-1)
        j2 = id2(n0)
        ! write(*,*) i1,i2,i3
        call triangle(x2(j1),y2(j1),x2(j2),y2(j2),x2(j3),y2(j3),u2,v2)

        q = .false.
!        puv = exp(-suv**2*(dd/1.0)**2/uradius2)
!        puv = exp( - (abs(u1 - u2) + abs(v1 - v2))/2.0)
!        puv = max(1.0_rp - abs(u1 - u2),1.0_rp) * max(1.0_rp - abs(v1 - v2),1.0_rp)
!        puv = (1.0_rp - abs(u1 - u2)) * (1.0_rp - abs(v1 - v2)) * (1.0_rp - abs(u2 - v2))
        puv = (1.0_rp - abs(u1 - u2)) * (1.0_rp - abs(v1 - v2)) 
!        call angle(x2(j1),y2(j1),x2(j2),y2(j2),x2(j3),y2(j3),p1,p2)
!        pang = exp( - (p1 - asix) - (p2 - asix))
!        pang = abs(p1 - asix)/asix + abs(p2 - asix)/asix

!        puv = 1.0_rp - sqrt((u1 - u2)**2 + (v1 - v2)**2)
!        puv = exp( - sqrt((u1 - u2)**2 + (v1 - v2)**2))
!        puv = puv * exp( - (1.0 - v2) - (1.0 - u2))

!        puv = puv*pang

!        if( v2 < 0.1 .or. u2 < 0.1 ) puv = 0.0
!        puv = puv * exp( - (1.0 - v2) - (1.0 - u2))
!        puv = puv * (1.0 - ((1.0 - v2) + (1.0 - u2))/2.0)

        if( pmin * puv >= fp0 ) then
!        if( puv > 0.0 ) then

           n0 = n0 + 1
           pmin = pmin *puv
           id1(n0) = i3
           id2(n0) = j3
           !              write(*,*) id1(1:n),i,'->',id2(1:n),j,s0,suv,ss,n0
           if( n0 < nmatch .and. i3 < size(x1) ) then

              call quest(i3+1,q)
              if( .not. q ) then
                 n0 = n0 - 1
                 pmin = pmin / puv
              end if

           else
              q = .true.
              if( pmin >= fp0 .and. n0 >= fn0 ) then
                 dmin = dist(x1,y1,id1(1:n0),x2,y2,id2(1:n0))
                 write(fmt,"(a,i0,a,i0,a)") "(",n0,"i3,a,",n0,"i3,f10.3,f8.5)"
!                 write(*,*) ":",trim(fmt),":"
!                 if( pmin > 0.1 ) write(*,fmt) '+',id1(1:n0),'->',id2(1:n0),pmin,n0,suvs(3:n0)
                 if( pmin > 0.0 ) write(*,fmt) id1(1:n0),' -> ',id2(1:n0),pmin,dmin
                 fn0 = n0
                 deallocate(fid1,fid2)
                 allocate(fid1(fn0),fid2(fn0))
                 fid1 = id1(1:n0)
                 fid2 = id2(1:n0)
                 fp0 = pmin

                 if( nitems < maxitems ) then
                    nitems = nitems + 1
                    items(nitems)%p = pmin
                    items(nitems)%id1 = id1
                    items(nitems)%id2 = id2
                 else
                    ! locate the most worse ones and replace it
                    pmax = 1.0
                    imax = 0
                    do i = 1,size(items)
                       if( items(i)%p < pmax ) then
                          pmax = items(i)%p
                          imax = i
                       end if
                    end do
                    items(imax)%p = pmin
                    items(imax)%id1 = id1
                    items(imax)%id2 = id2                    

                 end if
                 
              end if
           end if
        end if

     end if

  end do

end subroutine quest


function locate(d,id)

  logical :: locate
  integer, intent(in) :: d
  integer, dimension(:), intent(in) :: id
  integer :: i
  
  locate = .false.
  do i = 1, size(id)
     if( id(i) == d ) then
        locate = .true.
        return
     end if
  end do

end function locate

subroutine triangle(x1,y1,x2,y2,x3,y3,u,v)

  real(rp), parameter :: eps = epsilon(0.0_rp)

  real(rp), intent(in) :: x1,y1,x2,y2,x3,y3
  real(rp), intent(out) :: u,v

  real(rp), dimension(3) :: d
  real(rp) :: x
  integer :: i,j

  ! compute sides of triangle
  d(1) = (x1 - x2)**2 + (y1 - y2)**2
  d(2) = (x1 - x3)**2 + (y1 - y3)**2
  d(3) = (x2 - x3)**2 + (y2 - y3)**2

  ! sorting
  do i = 2, 3
     do j = 3, i, -1
        if( d(j-1) > d(j) )then

           x = d(j-1)
           d(j-1) = d(j)
           d(j) = x

        endif
     enddo
  enddo

  ! compute parameters
  if( d(3) > eps )then
     u = sqrt(d(1)/d(3))
     v = sqrt(d(2)/d(3))
  else
     ! it would be better put u,v to definition domain
     u = 0.0
     v = 0.0
  endif

end subroutine triangle


function dist(x1,y1,id1,x2,y2,id2)

  real(rp) :: dist
  integer, dimension(:), intent(in) :: id1, id2
  real(rp), dimension(:), intent(in) :: x1,y1,x2,y2
  integer :: i,n,k
  real(rp) :: d1,d2,mean
  real(rp), dimension(:), allocatable :: c

  n = size(id1)
  allocate(c(n-1))
  k = 0

  do i = 2, n

     d1 = sqrt((x1(id1(i-1)) - x1(id1(i)))**2 + (y1(id1(i-1)) - y1(id1(i)))**2)
     d2 = sqrt((x2(id2(i-1)) - x2(id2(i)))**2 + (y2(id2(i-1)) - y2(id2(i)))**2)

     if( d1 > 0.0_rp .and. d2 > 0.0_rp ) then
        c(i-1) = d2/d1
        k = k + 1
     end if

  end do

  if( k > 1 ) then
     mean = sum(c(1:k))/k
     dist = sqrt(sum((c(1:k) - mean)**2)/(k - 1.0))
  else
     dist = huge(dist)
  end if

  deallocate(c)

end function dist



function angels(x1,y1,id1,x2,y2,id2)

  real(rp) :: angels
  real(rp), parameter :: da = 6.28
  integer, dimension(:), intent(in) :: id1, id2
  real(rp), dimension(:), intent(in) :: x1,y1,x2,y2
  integer :: i,n
  real(rp) :: a1,a2,uu,vv
  real(rp), dimension(2) :: u,v
  
  n = 0
  angels = 0.0_rp
  do i = 2, size(id1) - 1
     u = (/ x1(id1(i)) - x1(id1(i-1)), y1(id1(i)) - y1(id1(i-1)) /)
     v = (/ x1(id1(i+1)) - x1(id1(i)), y1(id1(i+1)) - y1(id1(i)) /)
     uu = sqrt(sum(u**2))
     vv = sqrt(sum(v**2))
     if( uu > 0.0_rp .and. vv > 0.0_rp ) then
        u = u / uu
        v = v / vv
        a1 = acos(min(sum(u*v),1.0_rp))
        u = (/ x2(id2(i)) - x2(id2(i-1)), y2(id2(i)) - y2(id2(i-1)) /)
        v = (/ x2(id2(i+1)) - x2(id2(i)), y2(id2(i+1)) - y2(id2(i)) /)
        uu = sqrt(sum(u**2))
        vv = sqrt(sum(v**2))
        if( uu > 0.0_rp .and. vv > 0.0_rp ) then
           u = u / uu
           v = v / vv
           a2 = acos(min(sum(u*v),1.0_rp))
           angels = angels + abs(a1 - a2)
           n = n  + 1
        else
           angels = angels + da
        end if
     else
        angels = angels + da
     end if
  end do
  if( n > 0 ) then
     angels = angels / n
  else
     angels = da
  end if

end function angels


end module matcher

